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We study the effect of an artificial gauge field on the zero temperature phase diagram of extended 
Bose Hubbard model, that describes ultra cold atoms in optical lattices with long range interaction 
using strong coupling perturbation theory . We determine analytically the effect of the artificial 
gauge field on the density wave - supersolid (DW-SS) and the the Mott insulator-superfluid (MI 
-SF) transition boundary . The momentum distribution at these two transition boundaries is also 
calculated in this approach. It is shown that such momentum distribution which can be observed 
in time of flight measurement, reveals the symmetry of the gauge potential through the formation 
of magnetic Brillouin zone and clearly distinguishes between the DW-SS and MI-SF boundary. We 
also point out that in symmetric gauge the momentum distribution structure at these transition 
boundaries bears distinctive signatures of vortices in supersolid and superfluid phases. 
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I. INTRODUCTION 

The extended Bose Hubbard Model (eBHM) incorpo- 
rates the effect of long-range interaction upto various 
orders by adding the interaction terms between bosons 
localized at different lattice sites and has been widely 
studied [H-Q . Recently a number of ultra cold atomic sys- 
tems has been suggested as possible candidates to realize 
such models. This includes condensate of ultracold dipo- 
lar atoms like 52 Cr [l(J, heteronuclear polar molecules 
11] and Rydberg excited ultra cold atomic condensate 
12} . The effect of long range interaction can be mini- 
mally taken care of by adding a nearest neighbour inter- 
actions (NNI) in addition to onsite interactions to the 
prototype Bose Hubbard (BH) model The interest 
in this model stems from the appearance of new phases, 
namely the Density Wave (DW) phase and the Super- 
solid (SS) phase, apart from the Mott Insulator (MI) 
and Superfluid (SF) phase. Both DW and MI phases 
lack phase coherence as the superfluid order parameter 
vanishes and, both these phases are incompressible with 
a finite gap in the particle-hole excitation spectrum. 

For a system described by eBHM, in the limit where 
the hopping between neighbouring sites can be neglected 
in comparison to the onsite interaction, the system ad- 
mits DW and MI phases alternatingly as the chemical po- 
tential increases. As the hopping amplitude is increased 
in comparison to the onsite interaction and NNI the sys- 
tems makes a transition from DW to a SS phase or from 
MI to a SF phase depending on the chemical potential. 
The DW phase has alternating particle number at each 
site in contrast to the MI phase which has same number 
of particles at each site. This characteristic of DW phase 
is interesting to explore as this checkerboard arrange- 
ment continues even when the interactions are reduced, 
resulting in transition to SS phase prior to a SF phase. 

In the SS phase [T3 - [l6j the SF and crystalline orders 
co-exist resulting in the spatial modulation of superfluid 
density. This SS phase has been first claimed to be ex- 



perimentally observed in solid Helium (l7j . however the 
interpretation of the experimental results is a matter of 
continuing debate (18h21| |. On the other hand, the cold 
atomic condensates with long range interactions loaded 
in an optical lattice can act as a relatively reliable way 
to confirm the existence of SS phase. Particularly the 
effect of artificial gauge field, that can be created either 
by means of rotation [22l |23| or by imprinting motion- 
dependent laser induced phase on the internal states [24| 
on such systems, can lead to the formation of vortex 
states in such supersolid which has distinctive features 
as compared to the formation of similar vortices in an 
ordinary superfluid both in terms of critical velocity of 
superflow [25| as well as profile [26} . 

One of the most common probing technique to detect 
the different phases in BH model or eBHM is to study 
their momentum distribution. The information about 
the momentum distribution can be directly extracted 
from time-of-flight (TOF) absorption imaging in the long 
time limit of the freely expanding atoms released from 
the trap [13, HH . The experimental observation of MI- 
SF transition with ultra cold atoms in optical lattice [29| 
was based on this idea. In this context in the current 
work we study the effect of artificial gauge field on insu- 
lating DW and MI phases and their respective transition 
to SS and SF phases using a strong coupling perturbation 
approach. We analytically calculate the modification of 
phase boundaries due to the gauge field, momentum dis- 
tribution and its gauge potential dependence. We also 
compare the momentum distribution of the vortex states 
carying definite quasi angular momentum at the the DW- 
SS and MI-SF boundary. 

The effect of gauge field on the MI-SF transition in 
BH model has been studied extensively in recent times 
both within mean field approximation I30M36II and also 
by going beyond mean field description [37H40l | . Study 
of frustrated Bose-Hubbard model in presence of stag- 
gered fluxes using exact numerical diagonalization was 
also carried out recently and the dependence of the time 
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of flight image on the gauge potential was discussed [4l[ . 
In comparison, the DW - SS transition in eBHM in pres- 
ence of a finite flux due to the g aug e field (2(| as well as 
in presence of a staggered flux [42j is still in very early 
stage and was carried out only in mean field framework. 

The strong coupling perturbation expansion which we 
adopt here to study the eBHM model treats the hopping 
as perturbation |43lT45j ( for an alternative way of doing 
strong coupling expansion see [46j]) and results from such 
approach matches quite well with results from quantum 
monte carlo simulation [47j . Particularly momentum dis- 
tribution which can be compared with the experimental 
result of time of flight (TOF) imaging [48|, |49[ can be 
calculated using strong coupling approach even in the 
presence of an artificial gauge field where methods like 
quantum monte carlo is rather difficult to implement. 

The remainder of the paper is organized as follows. 
In section II, we present the model Hamiltonian and 
the formalism of our calculations using the strong 
coupling perturbation theory. In section III, we present 
the analytical expressions for the phase boundaries 
between the incompressible (DW and MI) phases and 
compressible (SS and SF) phases, in presence of artificial 
gauge field. Section IV includes the calculation of 
momentum distribution in presence of artificial gauge 
field obtained by using wave function calculated by 
strong coupling perturbation theory expansion and 
also, the quasi angular momentum distribution of 
the states. A brief summary of the conclusions is 
presented in section V and the other details of the per- 
turbation calculations are provided in Appendix A and B 



II. EXTENDED BOSE HUBBARD MODEL IN 
PRESENCE OF MAGNETIC FIELD 

To describe the effect of nearest neighbour interaction 
on ultra cold atomic condensate loaded in a square opti- 
cal lattice in the presence of uniform artificial magnetic 
field in transverse direction, we introduce the following 
extended Bosc Hubbard Hamiltonian. 

H = -^t i $b j + ~^ni{ni-l)-(i^ni+V^n i n j 

i.j i i i,j 

(i) 

The first term gives us the nearest neighbour hopping 
where the hopping matrix elements are non zero only 
for nearest neighbours and is given by tij — te l ^ ij with 
(f>ij — f r " dr.A(r) and A(r) is the vector potential cor- 
responding to the artificial gauge field. Here, b\, bi and 
fii are the boson creation, annihilation and number op- 
erators respectively. Here, V is the strength of nearest 
neighbour interaction that minimally captures the effect 
of long range interaction, \x is the chemical potential. We 
have rescaled the Hamiltonian by U and thus, all param- 
eters are measured in units of U. We neglect the effect of 
an overall trap potential assuming that it is sufficiently 



shallow and is neutralized by the effect of centrifugal force 
particularly at the central region of the condensate. We 
are particularly interested in the limit when Vd < U/2, 
where d is the dimension of the system which is 2 in the 
present case. In this limit, the alternating sites of the 
lattice in the DW phase contains no and no — 1 particles 
and such a phase is called no — \ DW phase. Along the 
t = axis the system will form alternative sequence of 
n — \ DW phases followed by MI phases of no particles 
at each site. In the rest of the paper that the alternative 
sites of DW phase have population ua and ub and finally 
set ua = no and ns = no — 1 to obtain the corresponding 
results for n — \ DW phase. 



A. Formalism 

Within the frame work of the strong coupling per- 
turbative expansion we calculate the ground-state en- 
ergy Eow{nA,nB) and Emi{tio) of the DW phase with 
ua and hb bosons on alternating lattice sites, and 
of the MI phase with ua — no bosons on each lat- 
tice site, respectively. We also calculate the energies 
of the DW particle-hole excitations and MI particle- 
hole excitations (states with an extra particle or hole), 
E p ^(n A ,n B ), E h D %{n A ,n B ), and £^ 7 r (n ), E^{n ), 
respectively. The unperturbed system corresponds to the 
case (t — 0). In this limit, ground state energy as well 
as the particle hole excitation energy can be analytically 
determined. Using the Rayleigh-Schroedinger perturba- 
tive expansion the ground state energy of DW and MI 
phases as well the energy of particle hole like excitations 
over these ground state is calculated in various orders of 
scaled hopping parameter t. 

Both the DW and MI states are gapped since the en- 
ergy to create a single particle-hole excitations is finite. 
With increasing t this gap energy starts decreasing. At 
the critical hopping parameter t — t c the energy to create 
a particle-hole pair vanishes and DW phase becomes de- 
generate with its particle and hole excited state. This 
gives the value of t at which DW-SS transition takes 
place. Thus, the phase boundary between the Density 
Wave and the Supersolid phase is determined by : 

E DW (nA,n B ) = E p ^ h °\nA,n B ) (2) 

Similarly, the phase boundary between the Mott Insula- 
tor phase and the Superfluid phase is determined as : 

E MI (n ) = EZ /ho \n ) (3) 

These conditions determine the particle and hole 
branches of both the insulating lobes (DW and MI), 
giving us u par and \i ho1 as functions of t 7 V 7 nA,riB (for 
DW phase) or t, V, uq (for MI phase). 
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B. Wave functions at zeroth order in t 

In this section we shall define the ground state wave 
functions for the DW and MI phases after setting the 
scaled hopping amplitude t = for hamiltonian defined 
in ([I]). In this limit these wavefunctions are determined 
by the competition between interaction energies alone 
and can be found out exactly. The wave functions for 
the particle and hole excited states above these ground 
states will also be mentioned both for DW and MI phase 
and we shall particularly emphasize the degeneracy asso- 
ciated with such excited states in the t — > limit. 

For the DW state, we divide the lattice into sublat- 
tices A and B, where each site in sublattice A contains 
ha particles and each site in sublattice B contains Ub 
particles. For t = 0, DW wave function can be written 
as : 



I* 



(0) \ 
DW I 




(4) 



where M is the total number of lattice sites, and |0) is the 
state with no particle. Here b\ ■ refer to boson creation 
operator on A and B sublattices, respectively. Since we 
are interested in calculating the wave function as well 
as energy of such a state for finite t as a perturbative 
expansion in the parameter i, the wavefunction defined 
in (j4j is also the wavefunction at the zeroeth order of this 
perturbative expansion. 

Unlike the ground state wave function defined in ((4]), 
the wave functions for the DW states with an extra par- 
ticle or hole for t = is degenerate. This is because 
when an extra particle or hole is added to the system, 
it can go to any of the M lattice sites. However in the 
case of a DW state, the alternatiing sites belong to A 
and B sublattice and contain different number of parti- 
cles. In the present case ua — and iib = uq — 1. 
Therefore a particle state over the DW ground state will 
consist of one particle added to any of site in sublattice 
B, which contains less number of particle in the ground 
state. All such states as well their linear combination 
are degenerate for t = 0. Since the total number of sites 
in B sublattice is 4r, therefore the dimension of this de- 
generate subspace of the one particle excitation is also 
4£. For the single hole type of excitation over the DW 
ground state similarly the hole can be created in any site 
that belongs to sublattice A, containing higher number 
of particle in the ground state. Therefore the dimension 
of the degenerate subspace of single hole like excitation 
is also 4^. Because of this degeneracy, to find out the 
states and energies of particle and hole like excitation for 
finite t as perturbative expansion in t we need to use the 
degenerate perturbation theory. 

To use degenerate perturbation theory to find out the 
wave function as well as the energy for particle or hole 
like excited state, we need to diagonalize the perturbation 
hamiltonian. To do this we write H given in |T]) as 



H — Hq + Hp 



(5) 



The unperturbed part, 

H = - ^2hi(hi - 1) - n^fii + V^hthj (6) 

i i i,j 

and the perturbed part, which is the kinetic energy (hop- 
ping) term. 



H f 



i,3 



(7) 



Now if we diagonalize Hp in the degenerate subspace of 
either particle or hole like excitation over the DW ground 
state, we shall find that the degeneracy is only lifted when 
we include the second order ( next nearest neighbour) 
hopping processes. This is because, each site in subltticc 
A(B) has only the sites of sublattice B{A) as its neigh- 
bour. Thus all the matrix element related by the first 
order hopping is and we need to go upto the second 
order to lift the degeneracy. Here we briefly mention the 
methodology. 

To find out the particle and hole excited state which 
will break the degeneracy when the perturbation is in- 
cluded we write it as a linear superposition of the degen- 
erate basis states, namely 



I* 



par(0)> _ 
DW I — 



M/2 



f WB b^%) 



I* 



hol(0)\ 
DW I 



M/2 



n A 



DWA b l* (0) ) 

DWI 



(8) 
(9) 



The correct choice for fj and will be obtained by diag- 
onalizing the second order perturbation due to the hop- 
ping matrix and identifying the corresponding min- 
imum eigenvalue. Therefore f® WB will be the eigen- 
vector of ^itjit^i with the minimum eigenvalue (e 2 £ 2 ) 
such that V,,* tj.l.yf 1 /" " = eH^!f WB and f l DWA 
will be the eigenvector of ^jtijt-^ with the mini- 



mum eigenvalue (e f) such that V, , ',,'//' 



DW A 



2+2 fDWA 



The normalization condition also requires 



that \f° WB \ 2 = 1 and Effi 2 \fP WA \ 2 = L 

Similarly for the MI phase, the non degenerate ground 

state wave function for t — is given by: 



M 



I* 



(0) \ 
Mil 



(10) 



fc=i 



where the index k refers to all the lattice sites. Here also, 
for t = the single particle or hole like excited states 
over this non-degenerate ground state is degenerate and 
dimension of these degenerate subspace is M, the total 
number of lattice sites. Since all sites are equivalent, this 
degeneracy is lifted in the first order correction of the 
degenerate perturbation theory, namely when the nearest 
neighbour hopping is included, unlike the previous case 
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in DW phase, where to lift the degeneracy one needs to 
include the next nearest neighbour hopping. Here very 
briefly we define the details. 

The single particle and hole excited state wave func- 
tions that will break the degeneracy for finite t is written 
as 



por(0)x 
MI I 



I* 



hol{0)\ _ 
MI I — 



1 



M 

V f MI 

-T / Jk u k 

k=l 



Mil 



M 
k=l 



(11) 

(12) 



The correct choice of the fk will be obtained by diago- 
nalizing the first order perturbation due to hopping and 
locating the minimum eigenvalue state. Thus fjf 1 is the 
eigenvector of the hopping matrix t kk > with the minimum 
eigenvalue (et) 

It may be noted that the Hp is same as the Harper 
hamiltonian whose spectrum is given by the Hofstadter 
butterfly. Therefore finding the solution that corresponds 
to the minimal eigenvalue of the hopping matrix tjj is 
identical to finding the band minimum in the Hofstadter 
problem and to locate their corresponding eigenstates. 
In the following section we shall provide the analytical 
expression of the perturbatively calculated energy of the 



ground state and the particle-hole excitation at finite t 
and evaluate the phase boundary of the DW-SS transi- 
tion and the MI-SF transition from this result. 



III. ANALYTIC EXPRESSIONS FOR THE 
MODIFICATION OF PHASE BOUNDARIES 

Using many body version of the Rayleigh-Schroedinger 
perturbation theory the ground state energy of the insu- 
lating phases ( DW or MI) as well as the energies of 
the particle or hole like excitation can be expressed as a 
power series in the scaled hopping amplitude t 



E(t) = t a E^ + t x E^ + t 2 E^ + t 3 E^ + 



(13) 



where E^ is the energy in the limit t = 0, and tE 
t 2 E {2) and i 3 ^ 3) are the first order, second order and 
third order corrections to energy. In the following we 
present such calculation upto the third order in t and the 
other details associated with the calculation are relegated 
to the Appendix A 

The expression of the ground state energy of the DW 
wave state per site is given as 



(i) 



J 



rptns 

M 



^o-l) 2 



zV 



n (n 



2n n 



H- 



V(z + 1) 



1 



V(z + 1) 



t 2 + 0(t A ) 



(14) 



r 



In the above expression z — 2d is the co-ordination num- 
ber of given lattice site which is in the case of a square 
lattice 4. The first term is the ground state energy at 
t = and the second term is the second order correc- 
tion due to finite t. Both the first and third order cor- 



rection vanishes, We now calculate the energies of the 
states with an extra particle or hole using Eqs.(j8]) and 
@ in the framework of degenerate perturbation theory. 
For the particle excited DW states, 



J 



E p 



E%fr+t [(no-i)+zVno-ii] 



+ t 



n (n + l)e 2 



n (n + l)z 



{1-zV) ((z-l)V) (2-(z + l)V) (l + (z-2)V) 
(n 2 -l)(e 2 -z) , 2n 2 (e 2 -z)^ 



(2 - zV) 



((z-2)V) 



0{t) 



(15) 



Comparing this expression with Eq. [13] we find that the manner the energy of the hole excited DW states can be 
first and third order correction vanishes. In a similar found out. This turns out to be 
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E h DW = E^+t Q {-(l + zV)(n -l)+[i} 
+ t 2r_«o(n - l)e 2 , n\t 2 



{ni-iy 



(1-zV) 
n (n - l)z 
V(l + (z-2)V) 



((z-l)V) (2-(z + l)V) 
(n%-l)(e 2 -z) , 2n 2 ( £ 2 -z) 



V(2-zV) 
I 



{\ + {z-2)vy 



o(t 4 



(16) 



In a similar way, the expression for energy of the MI 
phase and a single particle or hole excitation over it as 



J 



perturbative expansion in the scaled hopping parameter 
t upto third order is given by 



rpins i 



1) - fj.no + zV^- 



zt 2 Mno + l) +Q{ti) 



(1-V) 



(17) 



E^ s j + t°[n + zVn - M ] 
-*Vo + l)e + i 2 [n (n + 1) 



2 ^ , 2(z-e 2 ) 



2e 2 



(1-27) (1-V) 



n (n + 2) 



+ i 3 [-n (n + l)n 

- n (n + l)(n + 1) 

- n (n + l)(n + 1) 



(z - 2e)e + 



(e 2 -3z + 3)s 



{l-Vf 



2{l-VY 



, (2e 2 -6z + 6)e 2e(z - e 2 ) 2e(e 2 - 3z + 3) 4(z - 2)e 4e(e 2 - 3;z + 3) 



(s-e>- 

4e(e 2 - 3z + 3) 
(i7-7)(E/-27) 



(1-7) 2 



n (n + l)(n + 2) 



(1-27) 2 (1-V) 

e(z — 1) ze 



(1 - V) 4(1 - V) 2 



(1-27) (1-10(1-27) 
+ 0(0 (18) 



, /r = £jS? + t°(-K-l)-^Vn + u) 
- t[n e] +t 2 [(n + l)n 



2e 2 



t 3 [-?i (n + l)(n + 1) 



(z-2e)e + 



(1-27) (1-7) 
(e 2 -3z + 3)e' 

(i-v) 2 



- (n + l)(no - 1) 



2(1 -7) J 



n (n + l)n 



(z-6 2 ) £ - 



(2e 2 -6z + 6)e 2e(z - e 2 ) 2e(e 2 - 3z + 3) 4(z - 2)e 4e(e 2 - 3z + 3) 



n (n + l)(n + 1) 



(1-7) 2 (1-27) 2 (1-7) (1-27) (l-7)(l-27) 

4e(e 2 -3z + 3) \e(z-l) ze 



(U-V)(U-2V) 



-n (n + l)(n - 1) 



(1-7) 4(l-7) 2 



0(t q 



(19) 



r 



It may be seen that unlike in the case of DW phase, 
for the MI phase the the perturbative correction to the 
energy for finite t appears in the first order term itself and 
the second and third order corrections that appear in Eq.( 
[T8]) and (fl9|) is comparatively much smaller in magnitude. 
However we still kept the calculation upto the same order 
so that the entire phase diagram which comprises of DW 
as well as MI phase is given within the same order of 
the perturbation theory. In the following sections, using 
from the above expressions of energy of the ground state 
and their corresponding particle hole excitation we shall 



determine the MI-SF and DW-SS phase boundary from 
the relations ([2]) and ([3]). 



A. Determination of the DW-SS and MI-SF 
boundary from the strong coupling expansion 

Fig. ([1]) represents the phase diagram of the first four 
insulating lobes for the eBHM in two dimension (square 
lattice) in presence of artificial magnetic field. It shows 
the increasing stability of the insulating phase i.e the DW 
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and MI phases grows in size as the strength of magnetic 
field is increased from zero to finite values. This is due 
to the localizing effect of magnetic field on the moving 
bosons, thus favouring the insulating phases to occupy a 
larger area in the phase diagram. 

The minimum eigen value e used in calculation of en- 
ergy expressions, involves the diagonalization of the hop- 
ping matrix, which is gauge dependent. The location of 
the minimum eigen value e depends on the choice of the 
gauge potential, while the eigen value itself is not. We 
consider a rational flux v = p/q and calculate the mini- 
mum eigen values of the hopping matrix tij following [Hoj . 
The shape of the insulating lobes for the Mott phase is 
different from the DW phase, as observed in Fig ([lj) in 
two dimensions. This is because, the Mott states with an 
extra particle or hole, have corrections to first order in t 
, while the DW states have corrections to second order 
in t. So, as t — > 0, the slope of the Mott state will be 
finite, while it will vanish for the DW lobes. Hence, the 
shapes of the two insulating lobes are always different. 
Moreover, the shapes of the insulating lobes depend on 
the dimensionality of the system and also, on the appli- 
cation of artificial magnetic field [33] ■ The mean field 
theories always give a concave shape to the MI and DW 
lobes as the dimensionality only enters as a prefactor 
in the expression for [i as a function of t, while the strong 
coupling expansion easily distinguishes the shape of insu- 
lating lobes in different dimensions and also, in presence 
of artificial magnetic field. The determination of these 
transition boundary using strong coupling perturbation 
theory is one of the main results in this paper, Since the 
transition at the boundary in this model belongs to the 
universality class of 2 + 1 dimensional XY model, the 
critical exponents can also be found out through an ex- 
pansion of the chemical potential. This was done in (45j 
and we shall not discuss this issue further. 

In the next section, we calculate the effect of artificial 
magnetic field on the momentum distribution of the in- 
sulating phases (DW and MI), with both gauge choices, 
Landau gauge and symmetric gauge. 



IV. CALCULATION OF MOMENTUM 
DISTRIBUTION 

The standard experimental way of probing the proper- 
ties of an ultra cold atomic condensate is through TOF 
absorption imaging of the freely expanding atoms re- 
leased from the trap [28| . Same method is used for prob- 
ing the condensate loaded in optical lattice as well. The 
quantity that is measured experimentally in such TOF 
expansion is the momentum distribution of this ultra cold 
atomic ensembles. In this section we shall first provide 
the method of calculation of the momentum distribution 
within the framework of strong coupling expansion. Sub- 
sequently we provide our results and their analysis. We 
shall particulary show that in the presence of an artificial 
magnetic field the TOF immage actually depends on the 
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FIG. 1: (Color Online) Phase diagram for the eBHM in pres- 
ence of artificial magnetic field for Vd/U=0.2. The effect of 
increasing magnetic field is to increase the stability of the DW 
and MI lobes 

means to produce effective magnetic field. As already 
mentioned in all these discussion we have not taken into 
account the existence of an overall trap potential assum- 
ing it to shallow. Thus the discussion we present can be 
related to the bulk region of trapped condensate where 
the effect of the trap potential can be neglected. 

The momentum distribution n(k) is defined as the 
Fourier transform of the one body density matrix, 
p(r,r) = (r)ip(r )) , [52]. Therefore, 

n(k)= J dr J rfrV(r,r)e' k -( r - r ') (20) 

where V^( r ) an d ip( r ) are the bosonic field operators at 
the position r, respectively, and k is the momentum. The 
expectation value is taken in the many-body ground state 
and the corresponding wave function is determined using 
strong coupling perturbation theory as a power series in 
scaled hopping amplitude t ( in the unit of U). 

To evaluate (|20l) we expand the field operators in the 
basis set of Wannier functions, such that 

where M is the total number of lattice sites, w(r — R;) is 
the Wannier function localized at site I with position Ri 
and bi is the boson annihilation operator. Consequently, 
the momentum distribution becomes 

™( k ) = ^E^kO^kO^^e^- 1 ^ (21) 



where w(k) — J dru>(r)e lkr is the fourier transform of 
the usual Wannier function w(r). The summation in- 
dices le {A, B} and I e {A, B} includes the entire lattice. 
The wave functions for the insulating phases are calcu- 
lated pcrturbatively upto second orer in t which is the 
first significant order. The higher order correction can 
be neglected since i « 1, Upto second order, the wave 
function for the insulating state (MI/DW) can then be 
written as 



\tpi 



Here, |*°J 



K) + l€) + l€)+0(* 3 ) (22) 
' (0) \ defined in © and $Xj§ and 



I* 



I* 



I* 



(i) 



(2). 



MI/DW I 

£ (23, 



E 



(m\Hp\m){m\Hp\^> 



(0)\ 
insi 



E 0m > E 0m 



(24) 



In the expression for the first order correction (f23|) . the 
matrix element of Hp (hopping matrix) is taken between 
the first-order intermediate (excited) state |m > and the 
zeroth-order state defined in ((4]) and (fT0|) . where as \m ) 
in the expression of the second order correction (|24|) 
is the second order intermediate (excited) state. Also, 
Earn = Eq — E m is the energy difference between the 
first order intermediate state |m) and zeroth order state 
l^inl > an d E Qm > = Eq — E m > is the energy difference 
between the second order intermediate state |m > and 
zeroth order state 
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FIG. 2: (Color Online) Possibilities for the first order inter- 
mediate (excited) state for the DW phase in a 2x2 lattice. 
The figure (a) corresponds to the ground state wave function 
in the limit t = 



and contribute to the first order corrections in the ex- 
pression (|23|) . The further details of these calculations 
are given in the Appendix B. 

Next, we calculate the second order correction to the 
wave function of the DW state (eq |24|) . For every first 
order intermediate state |m), there will be number of 
possible second order intermediate states. For example, 
for the just discussed \m) = |n^' + — l,n^\n^) 

(Fig|2I&)) which is connected to the ground state by sin- 
gle hooping, the corresponding possibilities for \m ) are 
shown in Fig. (J3](a)-(g)) which is seven in number. 



To elucidate the meaning of the above terms appear- 
ing at the various orders of the perturbation theory, we 
depict the first and second order hopping processes re- 
spectively in Fig. ([2]) and Fig. ([3]) for the simplest case 
of a 2 x 2 square lattice unit cell. In the diagram given 
in Fig. ([2]) (a), we depict the ground state wave function 
for the DW state for t = for this unit cell whose wave 
function is can be written as 



|xT/(0) v _ I (1) „(2) (3) (4), 
I* DW/ — \ n A > n B > n A ' n B I 



(25) 



Here the superscript (i) correspond to the location of 
a given lattice site. All the subsequent diagrams in 
the same figure (Figj2j6)-(z)) depict all the possible |m) 
state defined in (|23|) which is connected to the DW 
ground state in Fig. [2][a) by a single hopping. For 
example, the state depicted in Fig. (0) correspond to 



+ l,n 



(2) 



l' n (^)> n (B) 



The matrix element 



of Hp between this state and the ground state can be 
calculated as 



(m\H P \*W) _ -Vn B (n A + l)t 



12 



E, 



On i 



(-2 + (z + l)V) 



(e) „„+i 



(a) », 
4 P~ 



16 



(b) n. 1 n A (C)» s+ 1 n, (d)„,- 
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FIG. 3: (Color Online) Second order intermediate states for 
jm >= + I, n^' — l,n^ 3 ',ng') ( the first diagram) 
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In the same way, for each possible first order inter- 
mediate (excited) state |m) (Fig. (U&)-(j)), the second 
order excited states are calculated and thus by adding 
all possible combinations, the second order correction to 
wave function is determined (Eq. [24l) . The details of the 
calculation is provided in Appendix B. 

With the help of above expressions for perturbation 



corrections, the DW ground state wave function can be 
determined upto the second order in perturbation theory 
using the general expression (|22|) and subsequently nor- 
malized within the same order of perturbation theory. 
Substituting of this wavefunction in the expression (|2ip 
for momentum distribution yields th n(k) for the DW 
phase as 



2n - 1 
2 

+ 2n 
+ 0(t 3 



,{z-l)V 



2(z-l) 



:V 2 



e(k) 



2- (z+l)V 

H - 1) 

2[2 - (z + 1)V} 2 ^ (z - l)V ^ [2 - {z + l)V] 



{nl - 1) 



[e 2 (k) -2dt 2 



(26) 



Similarly, for the MI phase with no particles on each 
lattice site, the momentum distribution in presence of 
gauge field is 



2n {n Q + 1) 



1 - V 



e(k)+n (n + l)(2n + l) 



x(e 2 (k) - 2dt 2 ) 



3-2V 



o(t 3 ) 



(27) 



In either of the above expression of momentum distri- 
bution (|26j) and (|27j) . the dispersion e(fc) is the minimum 
eigenvalue of the artificial magnetic flux dependent hop- 
ping matrix T or tij multiplied by a pre-factor jj for 
DW phase and jj for the MI phase where M is the total 
number of sites along a given direction. The matrix T in 
the lattice site basis can be written 



T 



hi 
hi 



tl2 ti3 
■ ■ ■ Uj 



1M 



tM,M-l tUM 



(28) 



As described in (JJ), Uj — te l ^ ij is the gauge dependent 
hopping amplitude from site i to site j and non- vanishing 
only if i and j are nearest neighbours. Through <f>ij — 
/ r * dr.A(r) this hopping amplitude explicitly depend on 
the gauge potential with A(r) as the vector potential. In 
the next sub section, we provide the results for the effect 
of artificial magnetic field on the momentum distribution 
of the DW and MI phases. 



Effect of the presence of gauge field on the 
momentum distribution 



The effect of artificial gauge field on the on the MI-SF 
and DW-SS transition boundary have distinctive features 
which can be demonstrated by plotting the momentum 



distribution derived in (|26j) and (|27j) in the k x — k y plane 
at these transition boundaries. Since the matrix T in 
(|2"5|) explicitly depends on the gauge potential, the mo- 
mentum distribution reflects in itself the gauge potential 
structure. 

To demonstrate this we calculate the dispersion e(fc) 
for two types of artificial gauge fields on the system, the 
Landau gauge and the Symmetric gauge and evaluate the 
corresponding momentum distribution structure. The 
vector potential in Landau ( L ) and Symmetric ( s ) gauge 
are respectively given as 

A L {v) = B(0,x,0) = 2-Kvxy (29) 
A s (r) = B(-y,x,0) = Tru(xy-yx) (30) 

where we consider that the flux due to the corresponding 
gauge field through the unit shell v = p/q is rational. 

In the Landau gauge A L , following [50] we denote the 
co-ordinate of a site i on the square lattice by pair of 
integers {m, n} in the unit of the lattice spacing a. For 
Landau gauge potential therefore the phase of the hop- 
ping parameter <f>ij — is if the link along x direction ( 
i — > j = {m,n} — > {m + l,n}) and (fiij = 2imv if the 
link is along the y direction such ( i — > j = {m, n} — > 
{m,n + 1} ). In terms of these notations the eigenvalue 
equation Tip = eip on the lattice can be written as 



^Pn— 1,771 ~t~ 

— l-Ktvri 



2niv n / 



i'n, 



m—1 



(31) 



The lattice wavefunctions that appears in Eq. (|31[) can 
be obtained by operating ip n ^ m by magnetic translation 
operator. Such operators are given as Tr = exp(jrR ■ 

[P + m^( r )]) [&3 where R is the lattice translational 
vector and it is known the operators along the x and y 
axis do not commute since T a xT a yT~^T~^ = exp(27rw). 
In case of Landau gauge and for v = ^ the required 
commutator is given by [T ga ^,T a y] — 0. To ensure that 
the wave function remain single valued at a given lattice 
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point as an unit cell is traversed, the enlarged unit cell 
known as magnetic unit cell therefore has gxl sites, as 
compared to the unit cell in the absence of such magnetic 
field. Correspondingly the Brillouin Zone (BZ) is reduced 
with — 7r < k y < tt, — n/q < k x < t: / 'q and is called 
the magnetic Brillouin zone (MBZ). As we shall see the 
momentum distribution n(k) shows the formation of such 
MBZ. 

Since [T, p y ] = in the Landau gauge the wave func- 



tion lpn,m{k x ,ky) 



1 f£ y I f l J-, I h' x 'I I , 



i(k x ,k y ) where k x , y are 



the components of the Bloch wave vectors. Substituting 
this expression in the eigenvalue equation (|31[) one gets 
the following one dimensional eigenvalue equation 

-t[e lk * (j) n+1 +e~ tk * cj) ni +2cos{ky+2Trvn)<j) n } = e<p n (32) 

The eigenvalues e(k) now can be determined from the 
condition 



dot 



,— ik 

vi i e 
e lfex ivi2 

: 



Mi e- 1 ^ 

•k- Ma 



M, 



N a 







9 -i e 







(33) 



with M„ = 2cos(k y a + 27m0) — e(fc) and, where 
n = 0,1, ...,q— 1. The eigenvalue matrix is q x q di- 
mensional because of the q times enhancement of the 
magnetic unit cell. And its minimum eigenvalue e(k) is 
q -fold degenerate since e(k x ,k y ) — e(k x ,k y + where 
a = 0, • • • , (q — 1) . Each of this minimum will corre- 
spond to a peak in the momentum distribution which is 
calculated using (|2"o]l and (|2"T)) . This has been plotted in 
Fig (Ufa)) and (b) for MI and DW phase. 

Thus in the Landau gauge potential, in MI phase 
one gets q peaks in n(k) at (k x ,k y ) = (0, ^p), with 
n = 0, 1, ..q — 1. Accordingly Fig (JIJa)) shows the 
momentum distribution for q = 4 which has peaks 
at (0,0), (0, ±7r/2) and (0,±tt). For the DW phase, 
v = 1/q with q = 4, we again get peaks at n(k) at (0, 0), 
(0, ±7r/2) and (0, ±7r), but in addition, there also exists 
small peaks at the corners of the reduced BZ, as shown 
in fig (jU[b)). The appearance of extra peaks at the BZ 
corners distinguishes DW phase from the MI phase. 
Since the DW phase consists of two interpenetrating 
sublattices A and B, even in the absence of any applied 
magnetic field, the BZ structure of the DW phase has 
double periodicity whereas the usual MI state has a 
single periodicity. This results in appearance of small 
extra peaks in the momentum distribution of DW even 
in the absence of magnetic field [49]. Upon application 
of magnetic field, we observe small peaks in the DW 
momentum distribution at the BZ corners, which is 
again attributed to the reduced periodicity of the DW 
phase compared to MI phase. 

When one uses a symmetric gauge potential given in 
(1301) the commutating magnetic translation operators 




■KlA 



t *** 



t 



FIG. 4: (Color Online) Momentum distribution for the (a) 
MI and (b) DW phase in landau gauge potential for v = 
j, plotted in the reduced brillouin zone — tt < k x < tt and 
-tt/4 < k y < 7T/4 



are T^qax and T^qay- The corresponding discretized 
eigenvalue equation was discussed in detail in [26| . 
Following the preceeding discussion it can be shown that 
the magnetic unit shell should be 2q x 2q of the unit 
cell in the absence field. The factor 2 comes because 
the phase accumulated when one goes around an unit 
cell in the square lattice is iriv in presence of the gauge 
potential (f30|) . Accordingly the MBZ will be defined 
as {-f q < k x < f q ,-f q < k y < §- q ). This results in 
the formation of peaks in momentum distribution at 
(±7rn/g, ±7rn/g). The momentum distribution of the 
MI phase in presence of symmetric gauge potential, for 
rational flux v = 1/4 is shown in fig ([5][a)). The peak 
in the momentum distribution is observed at the center 
((0, 0)) and smaller peaks are seen at the corners in the 
range (-7r/4 < k x < 7r/4, -7r/4 < k y < 7r/4). Same 
calculation for the DW phase, in a symmetric gauge 
potential, shows the existence small peaks at the BZ 
(doubly reduced) corners, in addition to the peak at 
(0, 0), as shown in fig ((5fb) ) . 

The effect of application of symmetric gauge potential 
to the system is equivalent to rotating the system (2(| [33[ 
and thus, leads to formation of vortices in the system 
when the superfluid order parameter is finite. Since we 
perform our calculations at the phase boundary, the re- 
sults for the momentum distribution of the DW and MI 
phases in presence of symmetric gauge potential, pro- 
vides information about the structure of vortex in a su- 
persolid (at the DW-SS phase boundary ) and vortex in 
a superfluid (at the MI-SF phase boundary). The pre- 
ceeding discussion demonstrates that rotating supersolid 
reflects some extra peaks in the momentum distribution 
in addition to the peaks observed in a rotating ordinary 
superfluid. These extra peaks are clearly demonstrated 
in fig © and is one of the central result of this work. 
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(a) (b) 




FIG. 5: ( Color Online) Momentum distribution for the (a) MI 
and (b) DW phase in symmetric gauge potential for v=l/4, 
plotted in the range — 7r/4 < k x < n/4 and — n/4 < k y < n/4 



Since momentum distribution can be measured us- 




FIG. 6: (Color Online) Extra small peaks in the DW mo- 
mentum distribution n(k) symmetric gauge potential. The 
arrows show the location of the extra peaks at intermediate 
points 



ing TOF absorption imaging technique, this provides a 
way to experimentally distinguish between the super- 
solid phase and superfiuid phase by comparing the re- 
spective vortex profile. To analyze this issue further in 
the next section we compare the momentum distribution 
n(k) in symmetric gauge to the momentum distribution 
corresponding to the many body states having definite 
quasi-angular momentum in a symmetric gauge poten- 
tial. This quasi-angular momentum are analogues of the 
Bloch-momentum, but for a rotationally invariant system 
and has been discussed in detail [53j 

B. Quasi-angular momentum distribution 

In this section we re evaluate the momentum distribu- 
tion n(k) but for a fixed quasi angular momentum and 
compare this with n(k) demonstrated in Fig. [5j 

To simplify calculation we approximate the wannier 
functions as delta functions to calculate the correlation 
function with a fixed value of angular momentum given 
to the system. The field operators can now be expanded 
as 

V>(r) = V e l2 * ml/M 5{r ~ n)6j (34) 

where m is the quasi-angular momentum, M is the num- 
ber of lattice sites. Substituting this in the expression 
pj]) yields 

n(k) = -J-V e i2ffm( '- i ' )/A/ J*(r-r i )<J(r-iy) 

x <i> ins \b}b l >\Tj} ins >j k -( r >- T ^ (35) 

The summation indices le {A, B} and I e {^4, B} includes 
the entire lattice. 

As compared to the expresssion (f2"Tj) , in (1551 the sys- 
tem has a prefixed value of quasi angular momentum 
and therefore we can associate a vortex phase with the 
system. 

For DW and MI phase, the quasi angular momentum 
distribution is given by : 



ra.Diy(k) 



M 



i.i' 

+2n 



V(z-l) 2U-V(z + 1)_ 
{nt - 1) 



2V 2 {z-l) 2 2[2U -V(z + l)} 2 UV(z-l) U[2U - V{z + 1)]. 



t 2 + 0(t 6 )) 



Similar expression for n(k) with a given quasi angular 
momentum for the MI phase can also be evaluated with- 
ing strong coupling expansion. For a square lattice with 
four sites, Fig. [7] (a) shows the momentum distribution 



for m = for DW phase. The m = state has a peak at 
k = 27m, with n as an integer. However for higher quasi 
angular momentum state m=l, n{k) vanishes at k = 2im 
in Fig. [7] (b) . Again the DW-SS phase boundary can be 
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FIG. 7: (Color Online) Quasi angular momentum distribu- 
tion of the DW phase for (a) m=0 ,(b) m=l, plotted over a 
range — 4-7T < fc^, fc M < 4-7T 



separated from the MI-SF phase boundary by noting the 
appearance of small extra peaks in the former case. To 
demonstrate these small peaks clearly we plot the cross 
sectional plots for the DW phase in Fig. |5J 

We also plot the quasi momentum distribution for the 
state m = 1 for MI and DW phases. 



coupling perturbation theory. Using the same approach 
we calculated the momentum distribution at this phase 
boundary which can be verified experimentally using the 
time of flight imaging. The momentum distribution re- 
flects the symmetry of the gauge potential. The momen- 
tum distribution shows distinctive feature at the DW-SS 
phase boundary as compared to the MI-SF phase bound- 
ary. Particulaly evaluating the momentum distribution 
for states with definite quasi-angular momentum at such 
phase boundary in a symmetric gauge we clearly demon- 
strate how this can be used to probe the vortex profile of 
the SS phase at the DW-SS boundary due to the action 
of the gauge field. There have been significant progress 
in cold atom experiments in identifying the enigmatic 
supersolid phase in ultra cold atomic systems [5J-l56|. 
Our calculation will hopefully augment further study in 
the behavior of such SS phase in presence of an artificial 
gauge field. 

RS is supported by a fellowship given by CSIR and 
the work of SG is supported by a grant from Planning 
section. IIT Delhi. 
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FIG. 8: (Color Online) The cross sectional plots for m=0 
and m=l quasi angular momentum state for MI((a),(c)) and 
DW ((b), (d)) phases. The DW phase shows distinctive peaks 
compared to MI phase 



Quasi angular momentum is connected with the phase 
information and hence, the vorticity. The fact that a zero 
quasi angular momentum state can be distinguished from 
a non-zero one by looking at the corresponding n(fc), al- 
lows the experimenters to verify that vorticity has en- 
tered in the system through the TOF measurement. 



V. SUMMARY AND CONCLUSION 



To summarize we have determined the modification 
of the DW-SS and MI-SF phase boundary correspond- 
ing to EBHM at zero temperature due to the effect of 
an artificial gauge field within the framework of strong 
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Appendix A: Perturbation corrections to the ground 
state and particle hole excitation energy in the DW 
and MI phases in presence of artificial magnetic field 



The first order correction is given as 



1. Corrections to grounds state energy 



Here we provide the detail expressions of 
^',En,En,En that appears in the equation 
T3| for DW phase. The zeroth order term is given by 



P (0) _ /, T .(o) ,„ , T (o) 



DW/ 



_ M , n A {n A ~ 1) + n B {n B ~ 1) , r n A n B 



n A + n B 
M o ) 



(Al) 



E 



(i) _ 

DW — 



(0) rr 1^(0) \ _ 



> = o 



This is because (n A — ^\n A ) 



(n B - l\n B ) 



(A2) 



0, where 



i and j can be any site. For the same reason all odd 
order correction like E\ also vanishes. The second order 
correction, which is also the first non- vanishing correction 
is given by 



E 



(2) _ \{^ DW {m A ,m B )\H P \^ w {n A ,n B ))^ 

dw — 2 t 



n A {n B + 1) 



n B (n A + 1) 



(n A — n B - 1) + (zn B - zn A + l)V (n B — n A — 1) + (zn A — zn B + l)V 

I 



Mzr 



(A3) 



Substituting the above expressions in the right hand side state upto the third order as 
of Eq. (fl~3|) one gets the ground state energy of the DW 



J 



EDw(n A , n B ) n A (n A - 1) + n B (n B - 1) 



M 



^n A n B n A + n B 

ZV 7L U ^ 



n A {n B + 1) 



n B (n A + 1) 



(n A — riB — 1) + (zn B - zn A + 1)V (n B — n A — 1) + (zn A - zn B + 1)V 



zt 2 + 0(t 4 ) 



r 



(A4) 



Setting n A = no and n B = no — 1 in the above expression 2. Calculation of energy of DW state with an extra 
we get the expression (|T4"1) . particle or hole 

Explicitly written wave function of the DW phase with 
an extra particle is 



Also, we can get the corresponding corrections for the 
ground state in MI phase by substituting in the above 
expressions for the DW state n A = n B = no- The ground 
state energy in the MI Phase given in Eq. (fT?)) can be 
obatined from the expression (|A4[) by putting n A — n B — 
no- 



\^ DW ^ 



M/2 

£ 

J 2 ) 



A ' A ■■ 



B 



B 



where /, 



DW{B) 



is the eigenvector of £\ tjit^i with 



the minimum eigenvalue (t 2 t 2 )ff WB , such that 



^ . . rDW(B) 



e 2 f 2 jDW(B) 
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The zeroth order energy of the DW phase is 



ppar(O) _ T par(0) 



par(O) 
DW 



> 



E 



(o) 



+ n B + zV'riA - /i 



First order correction in energy 



E 



par(l) 
DW 



par(0) 
DW 



> 



(A5) 



(A6) 



Again, all odd order terms vanishes. The second order 
correction to the energy of DW phase with an extra par- 
ticle is given as 



rp P ar(2) _ I < ^Dw( m A, m B)\Hp\ty P D J \n A ,n B ) > 

^DW — / j i w _ p s 

(n A + \){n B + l)e 2 t 2 

{n B - n A ) + (zn A - zn B )V {n A ■ 
n B (n A + l)e 2 t 2 



DW 



n A (n B + l)e 2 i 



2,2 



n B — 1) + (zn B - zn A + 1)V 
n A (n B + 2)zt 2 



(n B — n A — 1) + (zn A - zn B + l)V (n A - n B — 2) + (zn B - zn A + 2)V 
n B (n A + l)(e 2 - z)t 2 2n A (n B + l)(e 2 - z)t 2 



(n B - n A - 1) + (zn A - zn B )V [n A 

I 



n B — 1) + {zn B — zn A + 2)V 



(A7) 



Thus, adding equations (|A5|) . (|A7[) gives us the energy In the same way, the energy of the DW phase with an 
of DW state with an extra particle upto the third order. extra hole can be calculated and is given as 



E D ° w (n A ,n B ) = E l ^ w (n A ,n B ) - (n A - 1) - zVn B + \i + 

n A n B e 2 t 2 



(n B - n A ) + (zn A - zn B )V 
n B {n A + l)e 2 t 2 



n A (n B + l)e 2 t 2 



(n A - n B 



1) + {zn B - zn A + l)V 
(n A - l)(n B + l)zt 



(us - n A — 1) + {zn A - zn B 
n B {n A + l)(e 2 - z)t 2 



1)V 



(n A - n B 
2n A {n B - 



- 2) + {zn B - 
l){e 2 -z)t 2 



zn A 



{n B - n A - 1) + {zn A 



zn B )V 

r 



(n A — 71b — 1) + {zn B — zn A + 2)V 



2)V_ 
+ 0^) 



(A8) 



Again the expressions ([151) and (|T6|) are obtained by sub- 
stituting n A — uq and n B = uq — 1. 



where the coefficient fj are to be determined from the 



lowest eigenvalue Y^f tjf fy = ztf^ 11 as explained ear- 
lier. The zeroth order correction to the energy of MI 
phase with an extra particle is now given as 



3. Calculation of energy of MI state with an extra 
particle or hole 



The MI phase with an extra particle has the following 
wave function, 



M 

i*S (0) > = £/f>M 2) --4 J Vi,...4 M) > 

3 = 1 



E 



■par(0) 



MI 



= < * 



par(0) 
MI 



par(0) 
MI 



> 



= mi + n o + zVtiq - 



(A9) 



The first order correction to the energy of this state is 



7 par{\) 
J MI 



- par(0) I 
MI 



-.par(0) 
■ MI 



> = 



-(n + l)ei (A10) 
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The second order energy correction for the extra particle 
MI phase is as calculated below 



E 



,par(2) 

Ml ~ 



E 



< 



^Z r i(mo)\H P \^ r i (0 \n ) 



> 



E 



(2) 
MI 

no(no + 1) 

n (n a + 2)z 
2(1 -V) 



(K-E° m ) 

2(z-e 2 ) 



2c 2 



(1 - 2V) (1 - V) 

(All 



, P «r(3) <^ (0) K)|gp|^(^o ) >< *5E(«»o)|irp|*£K*b) >< *EJ(*o)|ffp|*£J("o) > 

MI - / , / , 



(E? n ~E°)(E° k -EV m) 



<K7i i( \n )\H P \*ZT(n a )>Y: 

(e 2 -3z + 3)e 



<^ m T>oMp|*m>'o)> 



-n (n + 1) < n 



- 2e)e + 



n ("o + 1) < (no + 1) 



(*-e 2 )e- 



(2e 2 -6z + 6)e 2e(z - e 2 ) 2e(e 2 - 3z + 3) 4(z - 2)e 



n (n + l)(no + 1) 



4e(e 2 -3z + 3) , 



(1-V) 2 (1 - 2V") 2 (l-V) 
e{z — 1) ze 



(1 - 2V) 



{l-V){l-2V) 



r-n (n + l)(n + 2) 



(l-V) 4(1 -V) 2 



(A12) 



Adding Eqs. dM|, (|AT0| . (fATTj) and (jAl2l) gives the 
energy expression for MI phase with an extra particle 
given in (|18[) . This expression recovers the known result 
of onsitc Bose Hubbard model when V 



[43|. 



The energy for the Mott phase with an extra hole is 
calculated, in the same way as above and is given in eq 
HE 



plicity). We apply perturbation theory on |^^^) to cal- 
culate I'I'^hV) up to the desired order cq([2"2"j). 
Using the diagrams given in Fig. the 1st order cor- 
rection to the ground state wave function of DW state in 
for a unit cell in a square lattice (for all possible |m)'s) 
can be explicitly written as 



Appendix B: Perturbative correction to the wave 
function and momentum distribution 

1. Corrections to the wave function for the 
insulating DW phase 

We start with the definition of the ground state wave 
function for the DW phase eq (l2~51) ( 2x2 lattice, for sim- 
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UT; (i) \ _ V- ~Eh' t W H^VIO) 

\^DWI — /-^i P I 



where 



-y/n B (n A + 1) u , (i) (2) (3) (4), . , (l) 1 (2) (3) (4) 
= L*i2|^ +1,« B - L,n A ,n B ) +ti 4 \n A + L,n B ,n A ,n B -1) 

4-f l«W „(2) i „(3) , „(4)v , |„(1) „(2) „(3) , 1 _(4) 1S1 

+ i32|« J 4,« s — ±,n A + i,n B ) + t 3i \n A ,n B ,n A + l,n B — i)\ 

\Jn A {n B + l) r . (1) (2) 1 (3) (4)> , | (l) (2) 1 (3) 1 (4), 
F2i|n A -l,n fl + L,n A ,n B ) + t 2 3\n A ,n B + i,n A - L,n B ) 

+ t^nf - l,nf,nf,nf + 1) + t^\n { l\n%\nf - l,n { B ] + 1)] 



Ei = (tie — n A — 1) + (ztia — zn B + 1)V 
E-2 = {nA - n B - 1) + (zn B - zn A + 1)V 



(Bl) 



To calculate the second order correction using eq ([M|) 
following Fig. Q we consider |m) 



1, n 



(2) 



The contribution of this |m ) state to the second order 
correction to DW wave function is 



(m\H f 



E, 



Om 



-\fn B (n A + l)t\ 
~E X 



The corresponding possibilities for |m ) for \m] 



l,n B — 1, n A ', n ( B ) are represented in Fig. dSJc)-(i)). 



(B2) 



For illustration let us choose 
Then 



2, n A \ n^ B 



E, 



Om 



= E - Ei - E { 



01 



2(n A - n B - 2) + (2zn A - 2zn B + 4)7 



\^D ] w) ai 



y/n B (n A + l)(n A + 2)(n B - 1) 



(#i)(£bi) 



2,nf 



'-12 



Similarly, we calculate the contribution from all of the 
\m ) states corresponding to this particular \m) . Sum- 
ming over this contributions we the second order correc- 
tion to the wave function due to this \m) state which 
looks like 



(.-C'lJi-C'OlJ 

+ ^WK + l)(^ + 2) ^ i4|n (, ) + _ _ i} 



*i2*32|ni 1 1J + l,ng J - 2,n A i) + l,n£>> 





(£l)(E 02 ) 




(rM - 


1- l)y/n B (n B -+ 


■1) 




(£i)(£ 03 ) 








1) 


(.El 



n B ^n A {n A + 
(Si) (£05) 



Vn^B(nA + l)(n B + l) 1(1),, (2) , (3) i^M.-u 

H , . , . t 12 t i3 \n A +L,n B - l,n A - L,n B + L) 

l-C'iJi-C'oeJ 

I + 1 (1) 1 (2) (3) (4) . 

+ fjp mtp \ t i2t34 n^ +l,n fl - i,n A ,n B - L) 
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Follwoing the same procedure for all other \m) states 
defined in Fig. @ and adding all contributions, we get 
the total second order correction to wave- function for the 
DW phase. 



where 



2. Terms in momentum distribution calculation 

The moemntum distribution in (|2ip can be written as 
|VF(k)| 2 



n(k) 



M 



■Y,(^Dw\blW\^ DW )e-^ R '- R ^ 



id' 



= n(k)(°> + n(k)W +n(k)W 



(k) (0) = ffi^^{^|6t6H^)}e- ik -^- R «': 
1,1' 



M 



£(<*MVl*£k) 



l,l' 



+ (*^I^VI*^»e- lk - (Ri - R '' ) 



,(k)( 2 ) 



l^(k)| 5 



A/ 



,(2) v 
DWl 



l,l' 



ik.(R ; -R ,) 



(2) 



(B3) 



(B4) 



(B5) 



Here we keep terms only upto second order and 
terms ^%\bjb f \^ w ), {^X\b\b r \^ w ) and 



J 



{*%\b\b^*%) are neglected, 
function from eq ([2"5j) we get g 



Using the wave 



and 



(0) _ n,-,i, : 2 f n A + n B 



i(k)W = - 



7i(k) w = |W(k)| 
n B (n A + l) \W(k)\ 2 



E 1 M/2 



R 



V 2 

n A (n B + l) |H/(k)| 2 



£ 2 



£*«. 



„ik.(R,-R,, 



M/2 ^ " 
Li' 



nB{n A + 1) + 1) 



Si 



£2 



e(k) 



M/2 



2£ 2 



2£f 



£1 



S 2 



(B6) 



(B7) 



(n fl +n A + l)[£ W ( »e ik '( R '- R .»'] 
ij'j" 



| W ^|2r n s( n ^ + !) , + !) ^b("a + 1) ^(n^ + 1) 1 ,-,,,,2,,,, .,./,2> 



2£ 2 



2£ 2 



£1 



x {n B +n A + l)(e 2 (k) - 2cft 2 ) 



(B8) 
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